%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%
%%  FFTfin
%%
%%  Authors: Dae Han Kang, Cong Han Lim, Justin Ormont,
%%  Michael Starr, Seeun Umboh, Irena Wang
%%
%%  This is a simple demonstration of the FFT algorithm
%%  described here: http://code.google.com/p/fft-group-representation-theory/
%%  Calculates the DFT of a delta function on a signal of length N.
%%
%%      
%%  This function has two parameters, p and k which describe the size of the
%%  signal for which the FFT will be calculated ( N = p^k)
%%
%%  Usage: FFTfin(p,k) 
%%  where p is prime > 2, and 0 < k < ~15 (depending on p)
%%  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function FFTfin(p,k),
% Initialize delta function
n=realpow(p,k);
f=zeros(1,n);
f(2) = 1;
f=f(:);

% Call FFT function with display
g=FFT_final(p,k,f);
end

% Calculates the intermediate transformation
function g=FFT_final(p,k,f)
n=realpow(p,k);
F=f;

% For each intermediate space
for j=0:k-1
    F=FFT_inner(p,k,j,F);
end
% Normalization
g=F/realpow(p,k/2);
% Plot result
mesh(real(g))
end

function G=FFT_inner(p,k,j,F)
  % Plot intermediate functions
  mesh(real(F))
  sleep(10.5)
  
  % Initialize return matrix
G=zeros(realpow(p,k-j-1),realpow(p,j+1));
n=realpow(p,k);

% Calculate coefficients
for r=0:realpow(p,k-j-1)-1
    for s=0:realpow(p,j+1)-1
        for a=0:p-1
            x=realpow(p,k-j-1)*a+r;
            y=mod(s,realpow(p,j));
            G(r+1,s+1)=G(r+1,s+1)+exp(-2*pi*i*(1/n)*((n+1)/2)*(realpow(p,k-j-1)*a*s+(s-y)*x))...
                *F(x+1,y+1);
        end
    end
end
end







